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Mechanically nonlinear energy harvesters driven by broadband vibrations modeled as white noise 
are investigated. We derive an upper bound on output power versus load resistance and show 
that, subject to mild restrictions that we make precise, the upper-bound performance can be ob- 
tained by a linear harvester with appropriate stiffness. Despite this, nonlinear harvesters can have 
implementation-related advantages. Based on the Kramers equation, we numerically obtain the 
output power at weak coupling for a selection of phenomenological elastic potentials and discuss 
their merits. 
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Energy harvesting from motion is a means to power 
wireless sensor nodes in constructions, machinery and on 
the human body [1,2]. A vibration energy harvester con- 
tains a proof mass whose relative motion with respect 
to a frame drives a transducer that generates electrical 
power. Linear resonant devices are superior when driven 
by harmonic vibrations at their resonant frequency, but 
perform poorly for off-resonance conditions. As real vi- 
brations may display a rich spectral content, sometimes 
of broadband nature, there has been considerable inter- 
est in using nonlinear suspensions to shape the spectrum 
of the harvester's response to better suit the vibrations 
[3-12]. The wider spectral response of nonlinear devices 
is expected to be beneficial for broadband vibrations. 

The studies so far indicate some advantages of nonlin- 
caritics for broad-banded vibrations, but little is known 
about which conditions make a nonlinear harvester fa- 
vorable compared to a linear one. This is due to lack 
of adequate theory and due to the studies being con- 
cerned about specific experimental or numerical exam- 
ples of nonlinear harvesters that are compared to specific 
examples of linear harvesters that could have been cho- 
sen differently. Furthermore, several studies do not con- 
sider the role of electrical loading which is known to have 
a dramatic influence on the consequences of mechanical 
nonlincaritics for the output power [13]. 

The Gaussian white noise process widely used in 
physics and engineering [14-17] is also important in 
studying broadband energy harvesting [4, 5, 13, 18-23]. 
It can model vibrations with a flat spectrum in the finite 
frequency range that the mechanical system can respond 
effectively to. The system itself then provides a frequency 
cut-off making the mathematical idealization of infinite 
bandwidth that characterizes white noise meaningful. 

Here we investigate theoretically the behavior of me- 
chanically nonlinear energy harvesters driven by a white 
noise acceleration. We derive rigorous upper bounds on 
the output power for arbitrary elastic potential and show 
that subject to mild restrictions on the device parame- 
ters, it is possible to find a linear device that performs 




FIG. 1. Vibration energy harvester model. 



equally well as the upper bound. We give a compact 
expression for the output power that we use to numeri- 
cally investigate the weak coupling limit of harvesters for 
different quartic polynomial potentials taking electrical 
loading fully into account. 

An energy harvester model that isn't technology spe- 
cific is shown in Fig. 1. The corresponding state space 
equations with a linear electromechanical transducer and 
a nonlinear mechanical suspension can be written 

x = v (1) 
mi) = -U'(x) - Tq/C - bv + ma, (2) 
-Rq = RI = V = Tx/C + q/C, (3) 

where m is the proof mass, x its relative displacement, 
v its velocity, U the open-circuit internal energy, q the 
transducer-electrode charge, V the output voltage, I the 
current, b the damping coefficient, R the load resistance, 
C the clamped capacitance and T is the transduction fac- 
tor. The device-frame acceleration y = —a is Gaussian 
white noise with a two-sided spectral density S a . The 
equations can represent a piezoelectric or an electrostatic 
energy harvester. An electromagnetic harvester gives the 
same mathematical structure, but different physical in- 
terpretation. We use charge as the independent variable 
[24]. Using voltage instead is physically equivalent and 
also common, sec e.g. [25]. 

Ensemble averages with respect to the stationary dis- 
tribution generated by the process (1-3) will be denoted 
by (•••). The mean output power P = (V 2 )/R will be 
our main object of interest. A number of other expres- 
sions for P immediately follow by using stationarity, (1) 
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and (3). We will use some of these expressions without 
giving the derivation. 

From (3), q = 0(T). The second term on the right 
hand side of (2) is 0(T 2 ) and can then be dropped in the 
limit r — > 0. This is the weak coupling limit, which in 
the stationary state has the reduced probability density 

W°(x, v) = cM-bv 2 /mS a - 2bU(x)/m 2 S a )/Z° t (4) 

where Z® t is a normalization constant [26]. We denote 
expectations in this limit by (■ ■ -)o. 

The important observation that the mean input power 
is Pj n = mS a /2 was made in [20] where it was proved for 
linear harvesters. For our nonlinear system and T = 0, 
all power is dissipated in the damper, (4) implies the 



equipartition theorem, and p„ = b(v 



mS a /2. For 



general T, consider the input energy J t ^ mavdt over a 
time interval. When the actually continuously differ- 
entiable a is modeled as white noise, the appropriate 
stochastic representation of the energy is a Stratonovich 
integral mvoda [27]. We have mvoda = mvda-\-mS a dt/2 
where vda is an Ito integral and has zero expectation 
[26, 28]. The input-energy expectation is then mS a dt/2 
which yields the stated expression for P n . 

The observation means that rj = 2P/mS a is an effi- 
ciency that should be maximized, as opposed to narrow- 
band harvesting where power transfer is maximized. It 
also implies a power balance 



P = mS a /2~b(v 2 }. 



(5) 



For linear harvesters, rj — > 1 as k 2 Q m — > oo where 
k 2 = T 2 /KC < I is the transducer electromechanical 
coupling factor, K is the open-circuit stiffness and Q m 
is the open-circuit quality factor[13]. Hence, it is impos- 
sible to improve significantly on a linear harvester that 
is already very efficient. The device in [29] for example, 
has k 2 Q m ss 7.8 resulting in rj rj 0.79. The great number 
of harvesters, especially those with small volume, that 
perform substantially below their theoretical maximum 
[2], suggests that the weak coupling regime nevertheless 
has great practical relevance. 

The load resistance determines the electrical time scale 
r = RC distinguishing different regimes of operation. 
When r is the fastest scale, i.e. r — > 0, we have [13] 

P ~ T 2 (v 2 )t/C ~ r 2 (v 2 ) r/C = T 2 TmS a /2bC. (6) 

From (3), it is readily proved that P = T 2 t{v 2 )/C — 
T 3 (i 2 )/C < T 2 t(v 2 )/C. One can also show that (v 2 ) < 
(v 2 )q. Hence, both asymptotic relations in (6) are upper 
bounds on the output power. We note that the bounds 
are valid for any U that permits a stationary distribution 
and that the output power is otherwise independent of U 
when t — ^ 0. 

When the electrical time scale is the slowest in the 
system, i.e. when r — ¥ oo, we have [5, 13] 

P ~ T 2 ((x - {x)) 2 )/rC ~ T 2 ((x - (,) ) 2 )o/rC. (7) 



The leftmost asymptotic formula in (7) is also an upper 
bound. This is seen by using (3) to find 

P = r 2 ((x - {x)) 2 )/rC ((q - (q)) 2 )/rC (8) 

which gives the inequality when dropping the second 
term. The rightmost asymptotic formula in (7) need not 
be an upper bound as can be inferred already from linear 
theory. We note that (7), in contrast to (6), is strongly 
dependent on U as it is proportional to ((x — (x)) 2 ). 

The maximum power as a function of r must necessar- 
ily be found at an intermediate value of r between the 
small-T and large-T regimes. Since the output power is 
respectively insensitive and sensitive to the nature of U 
in these two regimes, the degree to which the maximum 
power can be improved by mechanical nonlinearities is 
an open question. 

We first address the issue by deriving improved power 
bounds and comparing to linear behavior. Define z = 
q — (q) — D(x — (x)) — Bv and find the values of the 
constants B and D that minimize {z 2 }. Eliminate co- 
variances between x and q using T(xq) + (q 2 ) = and 
use P = T(qv) /C and (8) to write the minimum value as 



rC ((q~(q)) 2 ) 

r 2 ((x-(x)) 2 ) 



C^P 2 



P- 



(9) 



Next, use this to eliminate the variance of q in (8) and 
rearrange to obtain P = (1 — (z 2 ) /tCP)P u \ < P u i where 



Pm 



r 2 r(v 2 )((x-(x}) 2 ) 

C ((X-(X)) 2 )+T 2 (V 2 )' 



(10) 



We see that (10) agrees with (6) and (7) in their respec- 
tive limits and is a tighter bound. 

The rms frequency L0 m — y/ (v 2 ) /((x — (x)) 2 ) can be 
used to eliminate the displacement variance in (10). Us- 
ing P = mS a /2 — b(v 2 ) < P u i we find a lower bound 
on (v 2 ) which we substitute back into the power balance 
equation to obtain P < P U 2 where the new bound is 



Pu 



mS a 



T 2 r/Cb 



2 1 + T 2 T/Cb- 



(11) 



P U 2 is manifestly less than Pj n and is asymptotically ap- 
proaching the exact result at the extreme limits of r. 

The most optimistic estimate of output power permit- 
ted by (11) is found for load resistances such that w m r = 
1 and is P u2 ,Opt = (mS a T 2 /2Cb)/(2uj m + T 2 /Cb). The 
exact output power of an optimally loaded linear har- 
vester is P L in,o P t = {mS a T 2 /2Cb)/(2Lu + b/m + T 2 /Cb) 
where loq is the open circuit resonance. The two can be 
made equal by choosing loq = u> m — b/2m. The stability 
criterion k 2 < 1 dictates that such a linear device is re- 
alizable if uj m > b/2m + \T\/y/mC, i.e. unless we have 
very high damping. Therefore nonlinear harvesters are 
not fundamentally better than linear ones when driven by 
noise. 
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We now consider how to directly calculate the out- 
put power. From (1) and (3) it follows that V — 
(r/C) J_ oo cxp( — (t — ti)/r)v(ti)dti. Inserting this ex- 
pression into P = T(v(t)V(t)}, we obtain 

p2 poo p2 

P=-y o e-^(v(t)v(0))dt=-K vv (lM, (12) 

i.e. that the output power is proportional to the Laplace 
transform K vv of the velocity autocorrelation function. 

In the weak coupling limit r — > 0, we can approximate 
K vv by its value K° v for r = to obtain the leading 
order. K® v can be found from the transition probability 
by solving the Fokker-Planck equation corresponding to 
(1) and (2) with T = 0, i.e. the Kramers equation [30]. 

Without pursuing it further, we remark that an alter- 
native method to calculate the output power, and there- 
fore also K® v , would be to find a stationary solution of 
the Fokker-Planck equation for the energy harvester [13] 
in the weak coupling limit and use P ~ T(qv)o/C or 
P~T(vV) . 

We determine K® v numerically from the Kramers 
equation by orthogonal function expansions and matrix 
continued fraction methods following [30, 31]. The spa- 
tial basis functions are ip n (x) = \/W(x)ir n (x), n — 
0,1,... where W(x) = cyLp{-2bU(x)/m 2 S a )/Z , Z is a 
normalization constant and n n (x) are orthonormal poly- 
nomials with W(x) as weight function. We express all 
spatial-basis matrix elements in terms of the recurrence 
coefficients for 7r„ which are determined by adapting the 
Lanczos method described in [32] to continuous polyno- 
mials. Dimcnsionlcss variables distinguished by asterisk 
subscripts and based on a characteristic length scale l s 
and frequency scale oj s are used, e.g. P* = P/ml 2 Wg, 
= T/ y/mujgC, S a * = Safins and t* = oj s t. 

We first consider the much studied symmetric quartic 
potential U = Ax 2 /2 + Bx A /4, choose l s such that B* = 
Bl 2 /muj 2 = 1 and w s such that 7, = b/mui s = 1/100. For 
small values of r*, the output power in Fig. 2 collapses 
as predicted by (6) onto the same asymptotic form for 
both potentials. For S a * = 1.0 • 10~ 4 the mass vibrates 
around a potential minimum, giving a performance for 
larger r* that differs between the two cases due to their 
different linearized stiffness. At 5 a * = 0.1, the quartic 
term in the potential determines the behavior. In the 
intermediate case of S a * = 1.0 ■ 10 -3 , the two potentials 
give comparable maximum power even though there is a 
considerable difference between them for large t*. 

The upper bound (11) deviates significantly from the 
actual performance for the bi-stable potential at small 
S'a* when the longest time scale is that of interwell 
transitions as given by Kramers' rate problem [5, 14]. 
This demonstrates the necessity of the more complicated 
treatment, such as the numerical method used here, in 
predicting maximum power as opposed to bounding it. 

Fig. 3 shows the output power versus the parame- 
ter A,, = A/muj 2 when the load is optimized for ev- 
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FIG. 2. Output power P* versus electrical time scale r* for 
mono- and bistable potentials at weak coupling, 7, = 0.01, 
Sa, = 10~ 4 , 10~ 3 ,0.1 (bottom to top) and B, - 1. Open 
squares: Numerical solution for A* = —0.5. Solid circles: 
Numerical solution for A* = 0.5. Solid lines: Correspond- 
ing upper bounds. Dotted line: Solution from linearization 
around potential minima. 
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FIG. 3. Maximum output power P* as a function of A* at 
weak coupling, 7* = 0.01, P* = 1 and S a * = 10~ 3 . Solid 
circles: Numerical solution. Dashed lines: Solution from lin- 
earization around potential minima. Thin solid line: upper 
bound. Inset shows corresponding optimal load given by r*. 



ery A*. The value of the optimal r* in the inset varies 
correspondingly. The maximum power is obtained for a 
negative value of A*, i.e with a bi-stable potential, like 
demonstrated for a fixed load in [4]. But, as the bound 
shows, more power can be obtained for a linear device 
with ujq = u m . Also linear devices with the same stiffness 
as the nonlinear devices have at their potential minima 
can give more power. Therefore the motivation for uti- 
lizing nonlinear stiffness in noise harvesting is rather one 
of necessity than one of advantage. Implementation con- 
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FIG. 4. Quartic potentials. Dotted line: £ = 0, hardening 
Duffing spring; Dashed-dotted line: £ = y/2/3, negative tan- 
gential stiffness arises; Dashed line: £ = 2v2/3, bi-stability 
arises; £ = 1 symmetric bistable potential. 

straints such as, e.g., package size and/or beam dimen- 
sioning may prohibit linear operation. In this respect, 
we can think of the quartic term of the potential as a 
model of proof mass confinement or beam stretching at 
large amplitudes. 

We now consider a suspension made of a stable clastic 
material without built-in stress, choose U(Q) = and 
require U'(0) = 0, U"(0) > and U(x) > Vie ^ 0. 
The lowest order nontrivial polynomial form can then be 
parametrized as 

U{x) = ]-Kx 2 + ^L 3 + ^x 4 (13) 
w 2 yf2l 4Z 2 v 

where |£| < 1 and I is a length scale, see Fig. 4. Wc 
choose lo s = yJ K/m and l s = I as characteristic scales. 
A linear system with stiffness constrained to the same 
value K as in (13) is used for comparison. 

Fig. 5 shows that the nonlinear devices with f ^ 
give an S^-range of better performance than their linear 
counterpart. This is the case even with r* = 1 which 
is optimal only for the linear device. The consistently 
lower power for £ = is due to the stiffening nature 
of the potential which shifts the spectrum to higher fre- 
quencies. The other potentials have a range of softening 
behavior causing a shift to lower frequencies and higher 
power. The velocity spectral densities in Fig. 6 demon- 
strate this shift for the bistable potential with £ = 1 
and the monostable potential with £ = 2y/2/3. Despite 
their differences, these two potentials give very similar 
performance and spectral shape. Note that comparing 
instead to a linear system with luq = 0.5w s , which is 
within the S a * = 1 • lCP 3 -responses in Fig. 6, would 
give P*/S a *Tl ~ 50 outperforming all cases in Fig. 5. 
Even though we only considered simple phenomenolog- 
ical potentials (13), the broadening and flattening of 
the spectrum and the initial better-than-lincar power- 
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FIG. 5. Output power P* relative to acceleration spectral 
density S a * versus S a * both for r* = 1 (lines) and for optimal 
r* at each point (markers). Thin solid line: Linear device. 
7» = 0.01. 
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FIG. 6. Velocity spectral density versus frequency at weak 
coupling for S a * = 1 • 10~ 4 , 5 • 10 -4 , 1 ■ 10 -:! , 5 • 10~ 3 . Solid 
lines: f = 1. Dashed lines: £ = 2y/2/3. 



characteristic replicate experiments on a device with an 
asymmetric monostable potential [11]. 

We have shown that when driven by white noise, har- 
vesters with nonlinear stiffness do not have the funda- 
mental performance advantage over linear ones that one 
could have expected from their wider spectrum. This 
followed for efficient devices from considerations on in- 
put power and for general coupling from power bounds. 
Numerical examples were given for weak coupling. The 
findings do not preclude advantages of nonlinear-stiffness 
harvesters subject to vibrations significantly different 
from wide band noise, e.g. off-resonance, sufficiently 
band-limited vibrations. Implementation constraints 
may render a nonlinear stiffness unavoidable or a desired 
value of linear stiffness unattainable. We demonstrated 
advantages when linear stiffness was constrained. 
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